Percolation in living neural networks 
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We study living neural networks by measuring the neurons' response to a global electrical stim- 
ulation. Neural connectivity is lowered by reducing the synaptic strength, chemically blocking 
neurotransmitter receptors. We use a graph-theoretic approach to show that the connectivity un- 
dergoes a percolation transition. This occurs as the giant component disintegrates, characterized 
by a power law with critical exponent fi ~ 0.65. /3 is independent of the balance between excitatory 
and inhibitory neurons and indicates that the degree distribution is gaussian rather than scale free. 

PACS numbers: 87.18.Sn, 87.19.La, 64.60.Ak 
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Representing complex structures as connected graphs 
yields a simplification that retains much of their function- 
ality. It is therefore natural that network connectivity 
emerges as the fundamental feature determining statis- 
tical properties, including the existence of power laws, 
clustering coefficients and the small world phenomenon 
[l|. While experimental access to man-made networks 
such as the WWW or e-mail is feasible Q, biological 
ones such as the genetic networks must be painstakingly 
revealed node by node [3| . The connectivity in living neu- 
ral networks is even more difficult to uncover j^], since 
connections are hard to identify [3, Q and typically differ 
from brain to brain and from culture to culture [7|,|8|. Un- 
raveling the neural wiring diagram even in small cultures, 
with r^ 10^ neurons and '^ 10^ connections, is presently 
not feasible, although some progress has been attained 
in the study of the link between neural connectivity and 
information coding [5|, |9j, [1^ . 

Neural cultures derived from rat hippocampus develop 
into networks that display bursts of activity, governed by 
the presence of both excitatory and inhibitory neurons 
0, llll ■ In this Letter we present a novel experimental ap- 
proach to quantify statistical properties of such networks, 
and study them in terms of percolation on random graphs 
\1% [l3| . We control the connectivity of the entire net- 
work, gradually reducing the synaptic strength by means 
of chemical application. The initially connected network 
progressively breaks down into smaller clusters until a 
fully disconnected network is reached. The weakening of 
the network is quantified and the distribution of sizes of 
connected components determined by analyzing the neu- 
rons' response to stimulations applied simultaneously to 
the entire network. Viewed inversely, as the network's 
connectivity increases, a percolation transition occurs at 
a critical synaptic strength with the emergence of a gi- 
ant component, which increases as a power law with an 
exponent P ~ 0.65. 

Experimental setup and procedure.- Experiments were 
performed on primary cultures of rat hippocampal neu- 



rons, plated on glass coverslips following the procedure 
described by Papa et al. 'l4l| (Fig. [Ija)). The cultures 
were used 14-20 days after plating. The neural network 
includes N ~ 2.5 x 10^ neurons. The neural culture is 
placed in a chamber that contains two parallel platinum 
wires fixed at the bottom and separated by 15 mm (Fig. 
[IJb)). The neurons are electrically stimulated by apply- 
ing a 20 msec bipolar pulse through the wires. The cur- 
rent is controlled and gradually increased between subse- 
quent pulses, while the corresponding voltage drop V is 
measured with an oscilloscope. The chamber is mounted 
on a Zeiss inverted microscope with a lOX objective, and 
neuronal activity is monitored using the fluorescent cal- 
cium indicator Fluo-4. Images were captured with a 
cooled charge-coupled device (CCD) camera at a rate 
of 5 frames/sec, and processed to record the fluorescence 
intensity F of a sample of the culture including n ~ 400 
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FIG. 1: (a) Fluorescence image of a small region of the neural 
culture. Bright spots are cell bodies. Neural connections, 
mostly dendrites, are visible, (b) Sketch of the experimental 
setup, (c) Activity plot of the neural response. Black lines 
indicate those neurons that respond to the excitation, (d) 
F{t) signal for 3 neurons at increasing voltages. Vertical lines 
show the excitation time, and arrows the responding neurons. 



individual neurons as a function of time (Fig. [Ijd)). The 
images and the fluorescent signal are further analyzed to 
reject glia cells [l^- Neural spiking activity is detected 
as a sharp increase of the fluorescence intensity [la] . 

The connectivity of the network was gradually 
weakened by adding increasing amounts of CNQX 
(6-cyano-7-nitroquinoxaline-2,3-dione), the antagonist 
of the AMPA (alpha-amino-3-hydroxy-5-methyl-4- 
isoxazolepropionic acid) type receptors in the glutamate 
synapses of excitatory neurons. NMDA (N-methyl-D- 
aspartate) receptors were completely blocked with 20 
/iM of the corresponding antagonist APV (2-amino- 
5-phosphonovalerate) , enabling the study of network 
breakdown due solely to CNQX. To study the role of 
inhibition, we performed experiments with inhibitory 
neurons either active or blocked with 40 /iM of the 
GABA (Gamma-aminobutyric acid) receptor antagonist 
bicuculine. From here on, we label the network contain- 
ing both excitatory and inhibitory neurons by Gei, and 
the network with excitatory neurons only by Ge- 

The response of the network for a given CNQX concen- 
tration was measured as the fraction of neurons $ that 
fired in response to the electric stimulation at voltage V 
(Fig. Hfc)). Response curves $(t^) were obtained by in- 
creasing the stimulation voltage from 2 to 6 V in steps 
of 0.1 — 0.5 V. Between 6 and 10 response curves were 
measured per experiment, each at a different CNQX con- 
centration. Measurements were completed within 4 h, at 
the end of which the culture was washed of CNQX to ver- 
ify that the initial network connectivity was recovered. 

Model.- To elucidate the relation between the topol- 
ogy of the living neural network and the observed neural 
response, we consider a simplified model of the network 
in terms of bond-percolation on a graph. The neural net- 
work is represented by the directed graph G. Our main 
simplifying assumption is the following: A neuron has 
a probability / = f{y) to fire as a direct response to 
the externally applied electrical stimulus, and it always 
fires if any one of its input neurons fire. This ignores 
the fact that more than one input is needed to excite 
a neuron, and that connections are gradually weakened 
rather than abruptly removed. The model also ignores 
the presence of feedback loops and recurrent activity in 
the neural culture. However, we verified with numeri- 
cal simulations that relaxing these assumptions does not 
affect the validity of the model [1^ . 

Evidently, the firing probability $(/) increases with 
the connectivity of G, because any neuron along a di- 
rected path of inputs may fire and excite all the neu- 
rons downstream. All the upstream neurons that can 
thus excite a certain neuron define its input-cluster or 
excitation-basin. It is therefore convenient to express 
the firing probability as the sum over the probabilities ps 



of a neuron to have an input cluster of size s — 1 , 
$(/) = / + (1 — f)P (any input neuron fires) 

= f+{i-f)Y.Ps[i-{i-fr') 



(1) 



where we used the probability conservation X^sPs ~ 1- 
It is readily seen that $(/) increases monotonically with 
/ and ranges between <I>(0) = and $(1) = 1. The de- 
viation of <&(/) from linearity manifests the connectivity 
of the network (for disconnected neurons $(/) = /). Eq. 
(HI) indicates that the observed firing probability $(/) 
is actually one minus the generating function H(x) (or 
the z-transform) of the cluster-size probability Ps [l6|, 
H{x) — Y^^=i Ps^" — ^~^{f)i where x — 1 — /. One can 
extract from H{x) the input-cluster size probabilities p^, 
formally by the inverse z-transform, or more practically, 
in the experiment, by fitting H{x) to a polynomial in x. 
Once a giant component emerges the observed firing 
pattern is significantly altered. In an infinite network, the 
giant component always fires no matter what the firing 
probability / is. This is because even a very small / 
suffices to excite one of the infinitely many neurons that 
belong to the giant component. We account for this effect 
by splitting the neuron population into a fraction g that 
belongs to the giant component and always fires and the 
remaining fraction 1 — .g that belongs to finite clusters: 

*(/) = .9 + (1 - 5) [/ + (1 - f)P (any inp. neur. fires)] 

■DO 

= l-(l- g)Y^Ps{l-fr. (2) 

s=l 

As expected, at the limit of almost no self-excitation 
/ — )■ only the giant component fires, $(0) = g, and 
<&(/) monotonically increases to $(1) = 1. With a giant 
component present the relation between H(x) and the 
firing probability changes, obtaining 



iJ(a;) = ^p,x^ = l -$(/). 



(3) 



s=l 



In reality, the giant component is not infinite and it is 
measured from a sample which has n neurons. There- 
fore, it fires only after a non-zero, though small, firing 
probability /t is exceeded. To estimate this finite size 
threshold we note that when wc measure a giant compo- 
nent of size gn^ then the firing probability is 

^{f)^g[l^[l-fr),^g{l^e-f<^-). (4) 

This probability becomes significant at the threshold 

Ir^ignrK (5) 

Measured network response.- An example of the re- 
sponse curves ^{V) for a Gei network with n = 450 



neurons measured at 6 different concentrations of CNQX 
is shown in Fig. [S] At one extreme, with [CNQX] = 
the network is fully connected, and a few neurons with 
low firing threshold suffice to activate the entire cul- 
ture. This leads to a very sharp response curve that 
approaches a step function, where all neurons form a sin- 
gle cluster that comprises the entire network. At the 
other extreme, with high concentrations of CNQX (~ 
10 /LiM) the network is completely disconnected, the re- 
sponse curve rises moderately, and is given by the indi- 
vidual neurons' response. ^(V) for individual neurons 
(denoted as $oo(^)) is well described by an error func- 
tion <i>{V) = 0.5 + 0.5 erf (^^)- This indicates that 
the firing threshold of a neuron in the network follows a 
gaussian distribution with mean Vq and width 2ao. 
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FIG. 2: Response curves ^{V) for 6 concentrations of CNQX 
and n — 450 neurons. The grey bars show the size of the giant 
component. Lines are a guide to the eye except for 1 pM and 
10 /iM that are fits to error functions, with Vb = 3.96 V and 
CTo ~ 0.46 V for 10 /iM. Inset: Corresponding H{x) functions. 
The bar shows the size of the giant component for 300 nM. 



jump cannot be identified, or when the jump is compa- 
rable to the noise of the measurement, which is typically 
about 4% of the number of neurons measured. 

We studied the size of the giant component in a range 
of CNQX concentrations spanning almost three orders of 
magnitude from nM to 10 ^M in logarithmic scale. We 
define the control parameter c — 1/(1-1- [CNQX]/ii'c() as a 
measure of the synaptic strength, where the dissociation 
constant Kd = 300 nM is the concentration of CNQX at 
which 50% of the receptors are blocked TT'] . The param- 
eter c characterizes the connectivity probability between 
two neurons, and takes values between (full blocking, 
independent neurons) and 1 (full connectivity). Concep- 
tually, it quantifies the number of receptor molecules that 
are not bound by the antagonist CNQX and therefore are 
free to activate the synapse. 
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synaptic strength, c = 1/ (1 + [CNQX] / KJ 

FIG. 3: Size of the giant component as a function of the 
synaptic strength c, averaged over 18 Gei networks (•), and 
6 Ge networks (D). Lines are a guide to the eye. Some CNQX 
concentrations are indicated for clarity. Inset: Log-log plot of 
the power law fits <; ~ |1 — c/co| . The slope 0.65 corresponds 
to the average value of /3 for the two networks. 



Intermediate CNQX concentrations induce partial 
blocking of the synapses. As the network breaks up neu- 
rons receive on average fewer inputs and a stronger exci- 
tation has to be applied to light up the entire network. 
The response curves are gradually shifted to higher volt- 
ages as [CNQX] increases. Initially, some neurons break 
off into separated clusters, while a giant cluster still con- 
tains most of the remaining neurons. The response curves 
are then characterized by a big jump that corresponds 
to the biggest cluster {giant component), and two tails 
that correspond to smaller clusters of neurons with low 
or high firing threshold (Fig. [2]) . Error functions describe 
these tails well. Beyond these concentrations ([CNQX] 
> 500 nM for Gei networks) a giant component cannot 
be identified and the whole response curve is then also 
well described by an error function. 

Giant component. - The biggest cluster in the network 
characterizes the giant component, g. Experimentally, it 
is measured as the biggest fraction of neurons that fire 
together in response to the electric excitation. For each 
response curve, g is measured as the biggest jump A$, 
as shown by the grey bars in Fig. [2j The size of the giant 
component is considered to be zero when a characteristic 



Fig. |3] shows the breakdown of the network for both 
Gei and Ge networks. The giant component for Gei 
networks breaks down at much lower CNQX concentra- 
tions compared with Ge networks, and one can think of 
the effect of inhibition on the network as effectively re- 
ducing the number of inputs that a neuron receives on 
average. For Gei networks, the giant component van- 
ishes at [CNQX] ~ 600 nM, while for Ge networks the 
critical concentration is around 1000 nM. 

The behavior of the giant component indicates that 
the neural network undergoes a percolation transition, 
from a set of small, disconnected clusters of neurons to 
a giant cluster that contains most of the neurons. At 
the vicinity of the emergence of the giant component, 
the percolation transition is described by the power law 
g ^ \1 — c/co\^ ■ Power law fits for Gei and Ge networks 
give the same /? within the experimental error (inset of 
Fig. O, with Co = 0.36 ± 0.02, /3 = 0.66 ± 0.05 for Gei , 
and Co = 0.24 ± 0.02, /3 = 0.63 ± 0.05 for Ge ■ 

Cluster distribution analysis.- The construction of the 
experimental function H{x) defined in Eq. ([3]) allows the 
fit of a polynomial ^^ Psx'^ to determine the size distri- 
bution Ps{s) for clusters that do not belong to the giant 



component. Since / = $00 (V') is the response curve for 
individual neurons (Fig. ^ and x = 1 — /, the func- 
tion H{x) for each response curve is obtained by pfotting 
1 — ^{V) as a function of 1 — $oo(^)- For curves with 
a giant component present, its contribution is eliminated 
and the resulting curve normalized by the factor 1 — .g. 




FIG. 4: (a) Cluster size distribution Ps{s) for the experiment 
shown in Fig. [21 (b) Average cluster size, (s), and variance of 
the cluster size distribution, a^ = (s^) — (s)^, as a function of 
the concentration of CNQX, averaged over 15 Get networks, 
(c) Giant component firing threshold /t as a function of 1/g 
for two groups of experiments with 90 (D) and 400 (■) neu- 
rons measured, and for CNQX concentrations between and 
500 nM. Lines are least square fits. 



The inset of Fig. [2] shows H{x) for the response curves 
shown in the same figure. The corresponding Ps{s) dis- 
tributions, shown in Fig. HJa), are obtained from poly- 
nomial fits up to order 20 [15]. Since the cluster analy- 
sis is sensitive to experimental resolution, it is limited to 
g < 0.8, which corresponds to [CNQX] > 200 uM for Gei 
networks. Experimental resolution also limits the obser- 
vation of very small clusters for [CNQX] < 500 nM, since 
they are associated with values of ^(V) close to 1 [15| . 
Hence, the cluster size distribution of Fig.|31Ja) shows the 
correct behavior, but not the precise details. Overall, as 
shown in Fig. Hljb), the clusters start out relatively big 
and with a broad distribution in sizes, to rapidly become 
smaller with a narrow distribution for gradually higher 
concentrations of CNQX. Isolated peaks in Ps{s) indi- 
cate non tree-like clusters outside the giant component, 
in contrast to the model. This hints at the persistence 
of loops and at a strong local connectivity. While our 



sample covers only part of the culture, it does represent 
the statistics of the whole population. Sample size affects 
the noise level, but the overall cluster size distribution of 
Fig. m^a) is representative of the whole network. Our as- 
sumption that one input suffices to excite a neuron leads 
to an under-estimation of the cluster sizes [T^ , probably 
in direct proportion to the number of inputs needed to 
excite a neuron, which is on the order of ten [lO| . 

Finite size effects are observed in the behavior of the 
firing threshold for the giant component /t, which in- 
creases linearly with 1/g, as predicted in Eq. ([5]). frig) 
is measured for each concentration as the value of / at 
which the giant component fires. Fig. \^c) shows the re- 
sults for two groups of experiments with (n) = 90 and 
(n) = 400 neurons measured. Linear fits provide slopes 
~ 0.02 and ~ 0.005, which are of the same order of mag- 
nitude as 1/ri, with n the number of neurons measured. 

Discussion. - Neural cultures turn out to be an exper- 
imental systems in which a clear percolation transition 
can be mapped out in detail. The graph approach has 
proven remarkably successful in supplying a simplified 
picture for a highly complex neuronal culture, yielding 
quantitative measures of the connectivity. The measured 
exponent /3 appears to be independent of the culture de- 
tails, such as the ratio between excitation and inhibition 
or the variance between different cultures. Since /3 char- 
acterizes the distribution of connections per node, it is 
possible to estimate the connectivity distribution by com- 
paring the experimental value of /3 with the one obtained 
from simulations or theoretical developments. 

Numerical simulations of our model with a gaussian 
degree distribution provide /3 = 0.66 ± 0.05, as in the ex- 
periments 15]. Percolation on directed random graphs 
with power law degree distribution pfc(fc) ~ k~^ gives /3 
equal or larger than one, where its exact value depends 
on the degree exponent A [ig| . Since this value is clearly 
different from our experimental observations and simula- 
tions, we conclude that the connectivity distribution in 
the neural network is not a power law one. This may be a 
crucial difference between networks grown with cultured 
neurons versus those growing naturally in the brain. 
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